device,decomposed=0
loadct,5
!p.multi=0

nprocx = 8  ; phi
nprocy = 4  ; theta
nproc = 1
rds = 4.88e8 ; base of the convection zone in [m/s]
rr = 6.97e8  ; Solar radii in [m/s]
base = 'Netcdf_full_'

;;; READING COORDINATES
print,'reading model coordinates, r,theta and phi'
for j = 1, nprocy do begin  
    ext_j = string(j,format='(i02)') 
    for i = 1, nprocx do begin  
        ext_i = string(i,format='(i02)')
        ext_nproc = string(nproc,format='(i03)')   
        file = base+ext_i+'_'+ext_j+'_'+ext_nproc+'.nc'
        nproc++
        tmp = NCDF_OPEN(file)
        if (j EQ 1) then begin
            tmp = NCDF_OPEN(file)
            ncdf_varget,tmp, ncdf_varid(tmp, 'x'), phi_col
            if (i EQ 1) then begin
                phi = phi_col
            endif else begin
                phi = [phi,phi_col]
            endelse
        endif
    endfor
    ncdf_varget,tmp, ncdf_varid(tmp, 'y'), theta_row
    if (j EQ 1) then begin
        ncdf_varget,tmp, ncdf_varid(tmp, 'z'), z
        ncdf_varget,tmp, ncdf_varid(tmp, 'Time'), t
        theta = theta_row
    endif else begin
        theta = [theta,theta_row]
    endelse   
endfor
r = (rds + z)/rr  ;; fractional radius
theta = theta+!pi/2  ;; latitude

nr = (size(r))[1]
ntheta = (size(theta))[1]
nphi = (size(phi))[1]
print,'phi, theta, r', nphi,ntheta,nr
x = r#sin(theta)
y = r#cos(theta)

;;; READING VARIABLES
print,'reading model variables ...'
nproc = 1
for j = 1, nprocy do begin  
   ext_j = string(j,format='(i02)') 
   for i = 1, nprocx do begin
      ext_i = string(i,format='(i02)')   
      ext_nproc = string(nproc,format='(i03)')   
      file = base+ext_i+'_'+ext_j+'_'+ext_nproc+'.nc'
      print,'reading file: ',file
      nproc++
      tmp = NCDF_OPEN(file)
      ncdf_varget,tmp, ncdf_varid( tmp, 'u'), u_col
      ncdf_varget,tmp, ncdf_varid( tmp, 'v'), v_col
      ncdf_varget,tmp, ncdf_varid( tmp, 'w'), w_col
      ncdf_varget,tmp, ncdf_varid( tmp, 'oz'), oz_col
      ncdf_varget,tmp, ncdf_varid( tmp, 'Th'), Th_col
            ncdf_varget,tmp, ncdf_varid( tmp, 'p'), p_col
      if (i LE 1) then begin
         u_row = u_col
         v_row = v_col
         w_row = w_col
         oz_row = oz_col
         Th_row = Th_col
         p_row = p_col
      endif else begin
         u_row =  [u_row,u_col]
         v_row =  [v_row,v_col]
         w_row =  [w_row,w_col]
         oz_row = [oz_row,oz_col]
         Th_row =  [Th_row,Th_col]
         p_row =  [p_row,p_col]
      endelse
   endfor   
   if (j LE 1) then begin
      u = u_row
      v = v_row
      w = w_row
      oz = oz_row
      Th = Th_row
      p  = p_row
   endif else begin
      u = [[u],[u_row]]
      v = [[v],[v_row]]
      w = [[w],[w_row]]
      oz = [[oz],[oz_row]]
      Th = [[Th],[Th_row]]
      p  = [[p],[p_row]]

   endelse
endfor
print,'variables: '
print,'u',size(u)
print,'v',size(v)
print,'w',size(w)
print,'Th',size(Th)
print,'oz',size(oz)
print,'p',size(p)


end
